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Deep inelastic lepton-hadron scattering (DIS) has 
played a decisive role in our understanding of the deep 
structure of matter. The latest version of the experiment 
performed with colliding 27.5 GeV electrons or positrons 
and 820 GeV protons at HERA yields increasingly pre- 
cise data so that not only fully inclusive measurements 
can be used to study the physics of hadronic final states. 
In fact, the study of multi-jet events and event shapes 
has become an important project at HERA jj]. 

One of the important theoretical tools in the analysis 
of hadronic final states is perturbative Quantum Chro- 
modynamics (QCD). In order to make quantitative pre- 
dictions in perturbative QCD, it is essential to perform 
the computations (at least) at the next-to-leading order 
(NLO) accuracy. In the case of DIS such computations 
have so far been completed for one-jet inclusive and 2(+l 
beam)-jet cross sections O . In this letter we compute the 
3+1-jet exclusive cross section, defined using the k± al- 
gorithm ||, as a function of the jet resolution parameter 
at the NLO accuracy. We also compute various differ- 
ential distributions of the three-jet cross section, defined 
using the inclusive k± algorithm Q), that can be com- 
pared to recent data obtained by the HI collaboration 
M. With our Monte Carlo program one can compute 
the NLO cross section of any other infrared safe two- and 
three-jet quantity — the particular distributions shown 
here are given simply as illustration. 

An immediate application of our computation is, of 
course, the analysis of HERA multi-jet data. Using our 
program, it is also possible to find the higher order QCD 
correction to the forward jet cross section in the low-x 
regime in DIS, to which the O(a^) computation effec- 
tively gives only the leading-order prediction and large 
corrections may come from higher orders Our com- 
putation is also part of the 0(a^) two-jet cross section. 

In computing the NLO corrections we use the dipole 
formalism of Catani and Seymour Q| that we modify 
slightly in order to have a better control on the numerical 
computation. The main idea is to cut the phase space 
of the dipole subtraction terms as introduced in Ref . || , 
the details of applying it to the case of DIS will be given 
elsewhere. 



The advantages of using the dipole method arc the fol- 
lowings: i) no approximation is made; ii) the exact phase 
space factorization allows full control over the efficient 
generation of the phase space; iii) neither the use of color 
ordered subamplitudes, nor symmetrization, nor partial 
fractioning of the matrix elements is required; iv) Lorentz 
invariance is maintained, therefore, the switch between 
various frames (e.g. laboratory and Breit frames) can be 
achieved by simply transforming the momenta; v) the use 
of crossing functions is avoided; vi) it can be implemented 
in an actual program in a fully process independent way. 

In order to ensure the correctness of our results we 
checked the following points: (i) the correctness of all 
matrix elements and phase-space generation; (ii) the sub- 
traction term regularizes the real correction in all possible 
soft and collinear limits; (iii) the same quantity is sub- 
tracted from the real correction as added to the virtual 
one. 

We achieved full control on the matrix elements by 
utilizing the well-tested matrix elements that we had 
used in our program DEBRECEN || for computing NLO 
corrections to three- and four-jet production in electron- 
positron annihilation with different particles crossed 
into the initial state. In particular, we used the cross- 
ing symmetric tree-level five-parton helicity amplitudes 
of Ref. Q and the crossing symmetric three- and four- 
parton tree and one- loop amplitudes of Bern et al [^0) . 
We consider only virtual photon exchange, in which case 
the necessary one-loop four-parton plus a gauge boson 
matrix elements were also derived by Glover et al |pr| . 
Both sets of matrix elements were tested in e + e~ — > four 
jets NLO cross section computations and full agreement 
among the theoretical descriptions fl2]| and very good 
description of experimental data fis] ] were found. 

We have checked numerically that in all soft and 
collinear regions the difference of the real and subtrac- 
tion terms contain only integrable square-root singulari- 
ties. Furthermore, we have also checked that our results 
are independent of the parameter that controls the vol- 
ume of the cut dipole phase space, which ensures that 
indeed the same quantity has been subtracted from the 
real correction as added to the virtual one. 
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Finally, to have a further check of the computation, 
we utilized that the dipole method allows for the con- 
struction of a process independent programming of QCD 
jet cross sections at the NLO accuracy. We use the same 
program structure, with trivial modifications, to compute 
dijet and three-jet cross sections. Consequently, using 
the correct matrix elements and checking that our dijet 
cross sections are correct, we obtain the correctness of 
the three-jet result automatically. 

In order to check the strucure of our program we com- 
pared our NLO two-jet predictions to those of two other 
existing programs |T^,|l5|. In Ref. [p| a comparison of 
the these programs was presented using the NLO two- 
jet results for the modified JADE clustering scheme. In 
Table I. we recall the numbers presented there for one 
particular set of parton distribution functions together 
with the corresponding result of our computation. (We 
refer to Ref. jl5) for the precise meaning of these num- 
bers.) We find complete agreement with the predictions 
of DISASTER++. Our results also agree with the predic- 
tions of DISENT with slight tendency to lower values, but 
within statistical errors, from DISENT for small values of 
xb and the inelasticity y (first two bins). A more com- 
plete comparison will be presented elsewhere. 

TABLE I. Comparison of the NLO two-jet cross sections 
defined using the modified JADE algorithm and for MRSD'„ 
parton densities obtained with three partonic Monte Carlo 
programs DISENT, DISASTER++ and NL0JET++ (this work). 



bin 


DISENT 


DISATER++ 


NL0JET++ 


1 


578.4 ± 7.1 


585.0 ±2.6 


585.5 ± 2.0 


2 


361.1 ±3.5 


364.8 ± 1.5 


364.8 ± 1.5 


3 


120.1 ± 0.9 


119.1 ± 1.7 


122.9 ±0.5 


4 


95.4 ± 0.87 


98.1 ± 1.11 


97.7 ±0.6 


5 


54.9 ± 0.40 


55.3 ±0.46 


55.7 ±0.4 


6 


17.3 ±0.13 


17.5 ±0.06 


17.5 ±0.1 


7 


12.3 ± 0.15 


12.1 ±0.50 


12.4 ±0.07 


8 


8.52 ± 0.08 


8.61 ±0.12 


8.61 ±0.07 


9 


2.63 ± 0.02 


2.65 ±0.03 


2.63 ± 0.02 



In the case of DIS the observables may be defined ei- 
ther in the laboratory frame, in the hadronic center of 
mass frame (virtual boson plus proton rest frame), or in 
the Breit frame, which is characterized by purely space- 
like virtual boson momentum. In order to easily apply 
detector cuts, our program generates events in the labo- 
ratory frame. Due to the Lorentz invariance of the dipole 
method, we can freely transform the momenta to the 
other frames to compute the jet functions. Jet algorithms 
can be defined in any of these frames. The k± scheme we 
employ here is implemented also in the Breit frame ||. 
The full algorithm is a two-step procedure, the first one is 
a preclustcring of hadrons characterized by a pre-chosen 
stopping parameter c? C ut- For clustering the partons, we 
used the covariant _B-scheme. The result of the first step 



is the beam jet plus hard final-state jets. The second step 
is the resolution of the hard final-state jets into sub-jets 
characterized by a parameter y cut . We scaled the stop- 
ping parameter with the momentum transfer squared Q 2 
as d cu t = fcutQ 2 , and we chose y cu t — 1, when the second 
step is not carried out. Thus our three-jet cross section 
is a function of one parameter, / cu t- 

Once the phase-space integrations are carried out, we 
write the three-jet cross section (at a given / cut ) at NLO 
accuracy in the following form: 

ff3jct(/cut) = ^ / d-q f a (Ti, [i F ) a s (n R ) 2 B a (i], / cut ) 

a J 

+a s (n R ) 3 (fi( R ) (77, / cut , n R ) + B<p (ry, / cut , Mf ) 

+C a {r,,f mt ))]. (1) 

In Eq. ([j]) a s = a s /2ir, f a (f],jj,F) is the parton distri- 
bution function for parton type a at momentum frac- 
tion rj and factorization scale /ip. B a gives the Born 

( R) ( F) 

contribution, B a , B a , C a are the correction func- 
tions. In the higher order correction we separated the 
dependence on the renormalization scale /ir — xrQb.s. 
and factorization scale [If = xfQh.s., where Qh.s. is 
the hard scattering scale. Qh.s. is usually set event by 
event, therefore, the scale dependences can only be writ- 
ten in the factorized form if we consider differential cross 
section in Qh.s.- The B„ (?7, /cut, H>f) function is ob- 
tained by multiplying the Born squared matrix element 
with f} Q \nx 2 R , where (3q = | (HCa — 2A</), and inte- 
grating over the phase space. Apart from trivial factors 

( F) 

the B a (?7, /cut) A*f) function is obtained by convoluting 
the Altarelli-Parisi splitting functions with the integral of 
the colour-correlated Born squared matrix elements (see 
Eqs. (8.39) and (8.41) in Ref. §). The C a {r),f cut ) cor- 
rection function requires much more computer time to 
compute. 

Eq. (Q) shows that using the dipole method one may 
either compute the full cross section at the NLO accu- 
racy including the convolution with the parton distri- 
bution functions, or simply the parton level functions 

( F) 

B a {v, /cut), B a '(r), /cut) and C a (ri, /cut), which can then 
be convoluted with the parton densities after the Monte 
Carlo integration. The latter procedure is the proper one 
if we are interested in an a s measurement from DIS data 
(to avoid the reculculation of the Monte Carlo integrals 
after each step of the fitting iteration). However, it re- 
quires at least triple differential binning (one for / cut , one 
for the momentum fraction rj and one for the hard scat- 
tering scale Qh.s.) f°r the nine partonic functions (three 
for each of the three incoming flavor type), which is be- 
yond the scope of a letter. 

When presenting numerical results of three-jet cross 
sections in this paper we use the kinematic region used 
by the HI collaboration ||. For the basic DIS kinematic 
variables Q 2 , Xb and y = Q 2 /(sxb) we require 
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5 GeV 2 < Q 2 < 5000 GeV 2 , 

< x Bj < 1 , 0.2 < y < 0.6 . (2) 

Furthermore, we restrict the (pseudo)rapidity-range in 
the laboratory frame and the minimum transverse energy 
of the jets in the Breit frame as 

- K VietMb < 2.5 , E T >5GeV. (3) 

We choose the average transverse momentum of the jets, 

Q H .s. = ^E^(i) (4) 

as hard scattering scale. We also studied the other usual 
choice, when Q\ s is the negative of the momentum 
transfer carried by the virtual boson (Q 2 ), but have not 
found significant differences. 

In Fig. [|, we plotted the cross section convoluted with 
the CTEQ5M1 parton distribution functions Jl(| and us- 
ing the two-loop formula for the strong coupling, 

_ , > a s (M z ) ( ft \n(w(fi,M z )) 

a sW = TFT 1 ~ a s{Mz) 



w{n,M z ) 



0o w(n, M z ) 



(5) 



where w{q,q a ) = 1 - o a s (q o ) ln(q /q), a s (M z ) = 
0.118/27T and X = § [Y1C\ - 6C F T R N f - 10C A T R N f ), 
with Nf = 5 flavors. For the leading order results we 
used the CTEQ5L distributions and the one-loop a s 
(a s (M z ) = 0.127/27T and 0i = in Eq. @). We used 
running MS electromagnetic coupling at the scale of the 
virtual photon momentum squared. In order to show 
the weak dependence on the parton distribution func- 
tions, we plotted the perturbative predictions also using 
the MRST99 parton distribution functions [jl7| (dashed 
lines). There is no LO fit in the MRST99 set, therefore, 
in this latter case the LO curve is missing. 
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FIG. 1. The perturbative prediction for the three-jet cross 
section <T3j ct (/cut) as a function of / cu t. The solid lines were 
obtained using the CTEQ5 parton distribution functions and 
the dashed lines with the MRST99 set. 



The inset in Fig. [l] shows the K factor (ratio of the 
three-jet cross section at the NLO to that at the LO ac- 
curacy), indicating the relative size of the correction. The 
K factor exhibits a clear maximum, where the correction 
is more than 100 % and becomes smaller than one, imply- 
ing negative correction, for / cut < 0.2. The position and 
height of the maximum as well as that of the point K = 1 
depends on the lower limit on the Q 2 range. For instance, 
increasing the lower limit from 5 GeV 2 to 100 GeV 2 , the 
point K = 1 moves to about / cut = 0.03. Negative NLO 
correction indicates that the resummation of large loga- 
rithms is necessary in order to obtain a reliable prediction 
for smaller values of the stopping parameter. 

In Fig. H we study the scale dependences of the three- 
jet cross section at a fixed value of the resolution param- 
eter / cu t = 0.2. The strong dependence on the renor- 
malization scale observed at LO is significantly reduced. 
The factorization scale dependence is already not signif- 
icant at LO and does not change much. Setting the two 
scales equal, fi R — fip — A 4 , we can observe a wide plato 
peaking near fi — Qh.s.- 
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FIG. 2. The dependence of the three-jet cross section <73j e t 
at the fixed value of / cu t = 0.2 on the renormalization and 
factorization scales. 



The first measurement of the three-jet cross sections 
in DIS was published by the HI collaboration in Ref. ||. 
They used the inclusive k± algorithm to define the jets 
(the precise definition is given in Ref. j[8| ) , selected three- 
jet events and plotted differential distributions of certain 
kinematical variables. We computed the same distribu- 
tions at the LO and NLO accuracy and compared the 
predictions, corrected to hadron level (the correction fac- 
tors are about 0.8), to the published data. In Fig. || we 
present the dcr3j e t/dQ 2 differential distribution. We ob- 
serve that the LO prediction has a different shape than 
the data: too low for small values of Q 2 and too high at 
high values of Q 2 . The radiative corrections bring theory 
and experiment much closer: the NLO prediction gives 
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a remarkably good description of the corrected experi- 
mental data. Similar conclusions can be drawn from the 
distibutions of the other kinematical variables studied in 
Rcf. ||, but not shown here. 
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FIG. 3. The differential three-jet cross section for the 
5 < Q 2 < 5000 GeV 2 region compared to experimental 
data measured by the HI collaboration. The LO prediction 
(dashed line) was obtained with the CTEQ5L parton distri- 
bution functions and the NLO prediction (solid line) was ob- 
tained with the CTEQ5M1 parton distribution functions. 

In this letter we presented a NLO computation of the 
three-jet rate defined with the k± clustering algorithm in 
DIS. Our results were obtained using a partonic Monte 
Carlo program that is suitable for implementing any de- 
tector cuts. We found that the K factor is very sensitive 
to the allowed kinematic region. We demonstrated that 
the NLO corrections reduce the scale dependence signif- 
icantly. The NLO prediction is weakly dependent on the 
parton distribution functions. The NLO perturbative 
prediction gives a remarkably good description of data 
measured at HERA and corrected to parton level. The 
same program can be used for computing the QCD radia- 
tive corrections to the (differential) cross section of any 
kind of two-, or three-jet cross section or event-shape dis- 
tribution in DIS. We compared the two-jet rates obtained 
by our program to previous results and found agreement. 
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